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I. INTRODUCTION 

The variety of quantum Monte Carlo (QMC) meth- 
ods is the most universal tool for the numerical study of 
quantum many-body systems with strong correlations. 
So-called determinental Quantum Monte Carlo (QMC) 
scheme for fermionic systems appeared more than 20 
years agC''"-'^'''. This scheme has became standard for the 
numerical investigation of physical models with strong 
interactions, as well as for quantum chemistry and nano- 
electronics. Although the first numerical attempts were 
made for model Hamiltonians with local interaction, the 
real systems are described by the many-particle action 
of a general form. For example many non-local ma- 
trix elements of the Coulomb interaction do not vanish 
in the problems of quantum chemistry^ and solid state 
physics^. For realistic description of Kondo impurities 
like a cobalt atom on a metallic surface it is of crucial 
importance to use the spin and orbital rotationally in- 
variant Coulomb vertex in the non-perturbative inves- 
tigation of electronic structure. The recently developed 
dynamical mean- field theory (DMFT)ifor correlated ma- 
terials introduces a non-trivial frequency-dependent bath 
Green function. The extension^ of the theory deals with 
an interaction that is non-local in time. Moreover, the 
same frequency dependent single-electron Green-function 
and retarded electron-electron interaction naturally ap- 
pear in any electronic subsystem where the rest of system 
is integrated out. An interesting non-local effect due to 
off-diagonal exchange interactions may be responsible for 
the correlated superconductivity in the doped fuUerena^. 
It is worth noting that the exchange mechanism often 
has an indirect origin (like the super-exchange) and the 
exchange terms can therefore be retarded. 

The determinantal grand-canonical auxiliary-field 
schemei*Sii4 is extensively used for interacting fermions, 
since other known QMC schemes (like stochastic se- 
ries expansion in powers of Hamiltonian^^ or worm 
algorithms'^) suffer from an unacceptably bad sign prob- 
lem for this case. The following two points are essential 
for the determinantal QMC approach: first, the imag- 



inary time is artificially discretized, and the Hubbard- 
Stratonovich transformationi^ is performed to decou- 
ple the fermionic degrees of freedom. After the decou- 
pling, fermions can be integrated out, and Monte Carlo 
sampling should be performed in the space of auxiliary 
Hubbard-Stratonovich fields. Hirsch and Fye^ proposed 
a so-called discrete Hubbard-Stratonovich transforma- 
tion to improve the efficiency of original scheme. It is 
worth noting that for a system of N atoms the number 
of auxiliary fields scales cx N for the local (short-range) 
interaction and as iV^ for the long-range one. This makes 
the calculation rather ineffective for the non-local case. 
In fact the scheme is developed for the local interaction 
only. 

The problem of systematic error due to the time dis- 
cretization was addressed in several works. For bosonic 
quantum systems, the continuous time loop algorithmic, 
worm diagrammatic world line Monte Carlo schemei^, 
and continuous time path-integral QMGii overcame this 
problem. Recently a continuous-time modification of 
the fermionic QMC algorithm was proposecli^. It is 
based on a series expansion for the partition function 
in the powers of interaction. The scheme is free from 
time-discretization errors but the Hubbard-Stratonovich 
transformation is still invoked. Therefore the number of 
auxiliary fields scales similarly as the discrete scheme, so 
that the method remains local. 

The most serious problem of the QMC simulation 
for large systems and small temperatures is the sign 
problemi^ resulting in the exponential growth of com- 
putational time. This is a principal drawback of the 
QMC schemeii, but it is system dependent. For rela- 
tively small clusters, in particular for the local DMFT 
scheme, the sign problem is not crucialZii^. If we con- 
sider any subsystem obtained by integrating out the rest 
of the system, the Gaussian part as well as the interaction 
for the new effective action are non-local in space-time. 
Unfortunately, as we pointed out, the non- locality of the 
interaction hampers the calculation because it is hard to 
simulate systems with a large number of auxiliary spins. 
It is nearly impossible to simulate a system with interac- 
tions that are non-local also in time, when the number of 
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spins is proportional to (PN/St)"^ (/3 is inverse tempera- 
ture, and St is a time-slice). 

Recent developments in the field of interacting fermion 
systemai^ clearly require the construction of a new type 
of QMC scheme suitable for non-local, time-dependent 
interactions. In this paper we present a continuous-time 
quantum Monte Carlo (CT-QMC) algorithm which does 
not introduce any auxiliary-field variables. The princi- 
pal advantages of the present algorithm are related to 
the different scaling of the computational time for non- 
local interactions. The scheme is particularly suitable for 
general multi-orbital Coulomb interactions. The paper is 
aimed at a general description of the algorithm and the 
estimation of the computation complexity. We present 
the results for test systems to show an adequacy of the 
method. Moreover, an analysis of a non-trivial multi- 
band rotationally-invariant model with a time-dependent 
interaction is performed. This model demonstrates the 



main advantages of the numerical scheme. The paper is 
organized as follows. In Section II we discuss a general 
formalism. Section III contains several applications of 
CT-QMC scheme for simple systems in comparison with 
the exact solutions and results of the super-symmetric 
multi-band impurity problem. The conclusions are given 
in the Section IV. 



II. FORMALISM 

A. General principles 

One can consider the partition function for the system 
with a pair interaction in the most general case which 
has the following form: 



Z = TrTe-^, (1) 
S = 1 1 Kclddrdr' + J J J J WrfXcl,d^cl,d--'dridr[dr2dr'^. 



Here T is a time-ordering operator, r = {r, s,i} is a 
combination of the continuous imaginary-time variable 
T, spin orientation s and the discrete index i numbering 
the single-particle states in a lattice. Integration over dr 
implies the integral over r and the sum over all lattice 
states and spin projections: J dr = J2i J2s lo 

One can now split S into two parts: the unperturbed 
action 5*0 in a Gaussian form and an interaction W. We 



introduce as well an additional quantity a^.,, which can 
be in principle a function of time, spin, and the number 
of lattice state. The functions a^, will later help us to 
minimize the sign problem and to optimize the algorithm. 
Thus up to an additive constant we have: 



S = So + W, (2) 
So^JJ (< +// alliwlr^ +w;lr)dr2dr'2'^ cl,c^ drdr' , 

= / / / / ^nVl (4, c"-^ - <| )(c|;, c'-^ - <pdridridr2dr^. 

Now we switch to the interaction representation and make the perturbation series expansion for the partition 
function Z assuming Sq as an unperturbed action: 



oo oo 

^ = E ^fc= E I dri J dr[... J dr2k J dr'^^nk{ri,r[, ...,r2k,r'^k)^ (3) 

k=0 fc=0 

/t. ' 1' 2---' 2fc 



Here Zq = TrTe ^° is a partition function for the unper- Hereafter the triangle brackets denote the average over 
turbed system and the unperturbed system (for arbitrary operator A: < 

KKT =<T(ct,c''i-ap)-...-(c^, c'^--<-)>. (4) 
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Since 5*0 is Gaussian, one can 
apply Wick's theorem to transform iffl. Thus D^j ^^f' 

is a determinant of a 2fc x 2k matrix which consists of 
the two-point bare Green functions g^j., =< Tc\,c^ > at 
a^, = 0. Obviously, for non-zero a^, 

^"^■'■Z = ii^S;^ ~ J = 1, ••.2fc, (5) 
I 



where g^i{Ti, r'^, Tjj.) denotes the Green function for a 
general term of the series 

g;,(ri,ri,...,riJ = (D;,\-;;P)-ix (7) 

X < Tcld^icld-^ - ap) • ... • {cl ~ all') > . 

' ' 1 ' 1 '2k '2k 

Similarly, one can write formulas for other averages, for 
example the two-particle Green function. 

An important property of the above formulas is 
that the integrands stay unchanged under the permuta- 
tions ri,n>,n+i,ri>+i ^ rj,rj>,rj+i,rj>+i with any 
Therefore it is possible to introduce a quantity which 
we call "state of the system" and is a combination of the 
perturbation order k and an unnumbered set of k tetrades 
of coordinates. Now, denote ^Ik — klil.k, where the fac- 
tor fc! reflects all possible permutations of the arguments. 
For the Green functions, fc! in the nominator and denom- 
inator cancel each other, so that gK = gk- 

In this notation, 

Z^JnKD[K], (8) 
G;, = J gK^KDlK], 

where / D [K] means the summation by fc and integra- 
tion over all possible realizations of the above-mentioned 
unnumbered set at each fc. One can check that the facto- 
rial factors are indeed taken into account correctly with 
this definition. 



B. Convergence of the perturbation series 

It is important to notice that the series expansion for 
an exponent always converges for the finite fermionic 
systems. A mathematically rigorous proof can be con- 
structed for Hamiltonian systems. Indeed, the many- 
body fermionic Hamiltonians Hq and W have a finite 
number of eigenstates that is 2^, where N is the total 
number of electronic spin-orbitals in the system. Now one 
can observe that fifc < const ■ W^^^, where Wmax is the 
eigenvalue of W with a maximal modulus. This proves 
convergence of because the fc! in the denominator 



where 5ij is a delta-symbol. 



Now we can express the two-point Green function for 
the system using the perturbation series expansion 
(EJ. It reads: 



(6) 

I 

grows faster than the numerator. In our calculations for 
the non-Hamiltonian systems we also did not observe any 
indications of the divergence. 

The crucial point of the proof is the finiteness of the 
number of states in the system. This is a particular pecu- 
liarity of fermions. For bosons, on other hand, one deals 
with a Hilbert space of an infinite dimensionality. There- 
fore series like Q are known to be divergent even for the 
simplest case of a single classical anharmonic oscillatos^. 
It is important to keep this in mind for possible exten- 
sions of the algorithm to the electron-phonon system and 
to the field models, since these systems are characterized 
by an infinite-order phase space. 

It is also important to note that this convergence is 
related to the choice of the type of series expansion. In- 
deed, the series Q contains all diagrams, including dis- 
connected. In the analytical diagram-series expansion 
disconnected diagrams drop out of the calculation and 
the convergence radius for diagram-series expansion dif- 
fers from that of Eq.Q. 

For the purpose of real calculation, it is desirable to 
estimate which values of fc contribute the most to Z. It 
follows from the formula Q that 

< fc >=< W > . (9) 

This formula gives also a simple practical recipe for how 
to calculate <W >. For example, in an important case 
of the on-site Coulomb interaction, it gives information 
about the local density-density correlator. 



C. Random walk in if-space 

Although formula ^ looks rather formal, it ex- 
actly corresponds to the idea of the proposed CT-QMG 
scheme. We simulate a Markov random walk in a space 
of all possible K with a probability density Pk oc \^k\ 
to visit each state. If such a simulation is implemented, 
obviously 

g;, =i^/s (10) 



GL =Z-^< Tclde- 



E 

k 



j dri j dr[... J dr2fe5r'(^l: ^2fe)^fe(^l: ^2fc): 
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The overline here denotes a Monte C arlo averaging over 
the random walk, and s = riif/jfiif | is an average sign. 

Two kinds of trial steps are necessary: one should try 
either to increase or to decrease fc by 1, and, respectively, 
to add or to remove the corresponding tetrad of " coordi- 
nates" . 

Suppose that we perform incremental and decremental 
steps with an equal probability. Consider a detailed bal- 
ance between the states K and K' , where K' is obtained 
by an addition of certain tetrad r2k+i , ^2fe-i-i ' ''2fc-i-2 , f'2k+2 
to K. It should be noted that P(K) and P{K') appear 
under integrals of different dimensionality, respectively k 
and k + A. Therefore it is more correct to discuss the 
detailed balance between the state K and all K' with 
r2fe-i-i, r2fe_|_i, r2fc+2, '"2fc+2 corresponding to a certain do- 
main d^r. The detailed balance condition reads 

Pk-,k' Pk' d^r 



k-2 



k+2 



Pk'-, 



K 



Pk 



(11) 



where Pk^k' is a probability to arrive in K' after a 
single MC step from K. 

In the incremental steps the proposition for the 
four new points should be generated randomly. 
Denote the probability density in this generation 
p{r2k+i,r'2kj^i,r2k+2,r'2f.^2)- 1^ t^i^ ^l^^P is accepted with 
a conditional probability pk—^k' , then 



P> 



'k^K' ^ PK^K'Pir2k+i,r2k+i,r2k+2,r2k+2)d'^r- (12) 

For the decremental steps, it is natural to pick ran- 
domly one of the existing tetrades and consider its re- 
moval. So, 



'^K'^K = PK'^K/ik + 1). 



(13) 



Therefore, one obtains the condition for acceptance 
probabilities: 



Pk'^k 

PK-.K' 



K 



K' 



{k + l)p{r2k+i , r2k+i, r2k+2 , rafe+a)- 

(14) 



In principle, one can choose different piwrll'^.^rllXD j it 
is important only to preserve 114|l . We propose to use 



(15) 



IM^JJJjK'^'ldrdRdr'dR' 

to generate new points in the incremental steps. Then 
the standard Metropolis acceptance criterion can be con- 
structed using the ratio 



Ml 



ri...r2k+2 

,/ 

2k + 2 



(16) 



for the incremental steps and its inverse for the decre- 
mental ones. 

In general, one may want also to add-remove several 
tetrades simultaneously. A thus organized random walk 
is illustrated by Figure 1. The same Figure presents a 
typical distribution diagram for a perturbation order k 
in QMC calculation. 



Z — Zq -\- ... + Zk-2 + Zk-i + Zk + Zk+i -\- Zk+2 + 

k-^ k+^ 




FIG. 1: Schematic picture of random walks in the space of 
k; ri, ri, r2k,r'2k according to perturbation series expansion 
Q and an example of the histogram for the perturbation 
order k. 



D. A fast-update of Green function matrix 

The most time consuming operation of the algorithm is 
the calculation of the ratio of determinants and Green- 
function matrix. It's necessary for calculation of MC 
weights as well as for Green function. 

There exist so called fast-update formulas for calcu- 
lation of the ratio of determinants and Green-function 
matrix. Usual procedure takes operations, while the 
fast-update technique allows one to perform N"^ or less 
operations, where A^ is a matrix size. 

Our derivation of the fast-update formulas is a gener- 
alization of the Shermann- Morrison scheme for the deter- 
minatal QMC. Usually, the two types of steps (fc ^ fc-|-l 
and k — > fc — 1) are sufficient. However, the steps 
fc — > fc ± 2 can be also employed in certain cases (see later 
examples) , so we present here also formulas for that case. 

We use the following notation to derive the fast-update 
formulas: 

= Gi^nMnJ, (17) 

Hereafter summation over repeated indices is implied and 
(fc) denotes size of the matrix. In the last formulae matrix 
(fc) jg extended to be a (fc + 1) x (fc + 1) matrix with 
Mk+i,k+i = 1 and Mk+i,i = 0, Mi^k+\ = (it docs not 
change the ratio of determinants). Thus 

^(fc+i) ^ ^,^(fc)[i + AM^*^)]-!, (18) 
Using the standard 2x2 super-matrix manipulations one 
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can obtain the following expression for [1 + AM^'')] ^ matrix: 



[1 + am('=)]-i 



/ 1 + Gi,fc+iA 

G2,fc+iA~^-Rfc+i,i 



1 + G2,fe+iA~^-Rfe+i,2 

— X~^Rk+l,2 



Gi,fe+iA ^Rk+i.k 

G2,k+lX~^ Rk+l,k 



-G2,fc+iA~^ 



— A~^i?fe+i,fe A~^ 



(19) 



Then it's easy to the obtain fast-update formulas for the 
step A; + 1. Matrix M(*^+^) can be obtained from M^'^\ 



Finally the expressions for the matrix M*^*^"*"^^ and for the 
ratio of determinants have the following form: 



/ 



M' ■ 



— -^i,fc+iA ^ \ 



— -^fc,fe+iA ^ 
\~^Rk+i,k A~^ 



\ —A ^Rk+iA 

MIj = m1^j> +L,,k+iX-'Rk+i,j 



(20) 



detDC'+i) 



detD(fe) - Gfe+i,fe+i - Gk+l,iM^^)Gj^k+l - A - , 



where «, j = 1, A;. For the step k — 1 (removal of the 
column and row n) the fast update formulas for matrix 
j^{k-i) g^jjjj ^j^g ratio of determinants are as follows: 



M, 



(fc) _ 



^,3 



(21) 



dctM''°' 
detMC'-i) 



(fc) 



dctD^*-) 

One can also obtain fast-update formulas in the same 
I 



manner for steps fc ± 2. Let's introduce a 2 x 2 matrix A: 



G. 



1,1 



Gn iMi jG 



J ^],1' ■ 



(22) 



where q,q' = k + l^k + 2. Then the fast-update formulas 
for a step k + 2 look like 



j\4-(fe+2) = 



^^k+l.q'^q'A 

V ^K+2,q'^q'A 





~-^l.<if\,it+l 




\ 




— Lk,qXq^k+l 


— Lk,qXq^k+2 






\ -1 

^fc+l,fc+l 


^fe+l,fe+2 




~^k+2,q'^q' ,k 




'^fc-|-2,fe+2 


/ 



M' ■ = M.^*:^ 



^i,q\,q'^l\j^ 



detD(M 

I 



(23) 



where i,j = l,...,k. For the step k — 2 (removal of two 
columns and two rows n + l,n + 2) matrix A has the 
following form: 



where q,q' = n+l,n + 2. Then the fast update formulas 
for the matrix M^*^"^^ and the ratio of determinants are 



A, 



9i9' 



M„ 



(24) 
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as follows: 

Mg-') = - J A- J, M^fJ , (25) 

dctD"°~^' _ dctM<'°' _ HptTM 

Using the fast update formula for ratio of determi- 
nants, the Green function can be obtained both in imag- 
inary time and at Matsubara frequencies: 



9r' = ffjr' - E ffJr>^*,i5or" 



(26) 



2J 

Here goii^) is a bare Green function. 

Higher correlators can be obtained from Wick's theo- 
rem, just as in the auxiliary-field quantum Monte Carlo^ 
scheme. Also note that it's convenient to keep in mem- 
ory only the inverse matrices M instead of direct D in 
simulations. 



E. The sign problem 

A proper choice of a can completely suppress the sign 
problem in certain cases. To be concrete, let us consider 
a Hubbard model. In this model the interaction is local 
in time and space, and only electrons with opposite spins 
interact. Therefore it is reasonable to take — 5{t — 



T')6{i — i')a^, similar for ai, and a| 
perturbation W becomes 



a 



0. The 



,bba 



d = U {n^{T)-a^){n^{T)-ai)dt (27) 



Here the Hubbard U and the occupation number oper- 
ator n — c are introduced. The Gaussian part of the 
Hubbard action is spin-independent and does not rotate 

I T 

spins. This means that only g|, do not vanish, and the 
determinant in jSJl is factorized 



jjrir2 



rir3...r2k~l jjr2r4...r2k 



'^l^2---^2k 



^l^3---'''2k~l '2'4---'2fc 



r,r; ...ri 



(28) 



For the case of attraction [/ < one should choose 



(29) 



I T 

where a is a real number. For this choice = ffj , and 
therefore = Z^i. is always positive in this case, as 
follows from formula 

This choice of a is useless for a system with repul- 
sion, however. Compared to the case of attraction, an- 
other sign of w at Qfj = ai results in alternating signs of 
il/c with odd and even Another condition for a is 
required. The particle-hole symmetry can be exploited 
for the Hubbard model at half-filling. In this case, the 
transformation cj — s- C| converts the Hamiltonian with 
repulsion to the same but with attraction. Therefore the 
series Q in powers oiW = U ^ — Q!)(^i(''') — (y)dt 



does not contain negative numbers, in accordance to 
the previous paragraph. The value of the trace in © 
is independent of a particular representation. In the 
original (untransformed) basis the above W reads as 
U j[n^[T) — Q;)(n|(r) — 1 -|- ajdt. We conclude that 



1 



(30) 



eliminates the sign problem for repulsive systems with a 
particle-hole symmetry. Of course, the average sign for a 
system with repulsion is not equal to unity in a general 
case. 

It is useful to analyze a toy single-atom Hubbard model 
to get a feeling for the behavior of the series (|3J|. The two 
parts of the action are 

5o = /(-Ai + Ua{)n^{T) + + U a^)ni{T))dT-{?,l) 
W = U finiir) - aj.)(n|(r) - a|)dr. 

Here /i is a chemical potential. For a half-filled system 
/i = U/2. For this model, it is easy to calculate terms of 
the series for Z explicitly. We obtain 

X (^l + e'3(f-C'"T)(l-a^i)'=) . (32) 

Consider the case of repulsion {U > 0). Let us use the 
condition H3()|l for an arbitrary filling factor. The later 
expression can be presented in the form 

X (l + e/5(-A'+t^")(l (33) 

For /i = C//2 the value of ilk is positive for any a. For a 
general filling factor, the situation depends on the value 
of a. For < a < 1 negative numbers can occur at 
certain k. Outside this interval all terms are positive, 
and there is no sign problem for the single-atom system 
under consideration. 

Since the sign problem exists already for the impurity 
problem for < a < 1, such a choice is also not suitable 
for the iV-atom repulsive Hubbard system. On the other 
hand, minimization of W requires a to be as close to this 
interval as possible. Therefore it is reasonable to take 
a = 1 or slightly above. This is the same as zero or 
a small negative value, since — 1 — ai- We use the 
similar choice of a's for more complicated multi-orbital 
models and always obtain a reasonable average sign. 

Finally, one may prefer to have a perturbation that is 
symmetrical in spin projections. Formula H29|) for the 
attractive interaction is already symmetrical. For the 
case of repulsion we propose to use a symmetrized form 

^(n| + a){ni ~ 1 ~ a) + j{n^ - 1 - a)(n| + a) (34) 

with some small positive a. 

There is another argumentation why the presence of 
Qf's in Eq. (|34|l is very important. Indeed, proper choice 
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of a make the average of Eq. (|34|l negative. We can call 
such an interaction "virtually attractive in average". It 
makes possible to obtain the A;-series with the all-positive 
integrals in the expansion, whereas the same series with- 
out a's is useless due to the alternative signs of integrals. 
We believe that the similar reasoning is valid for the non- 
local interaction. Note however that the proper choice of 
the a's depends on the particular system under calcula- 
tion. For now, we cannot offer a general recipe. In a 
certain situation, the expressions under the integrals are 
not always positive, and the exponential falloff occurs 
for the large systems or small temperature. The practi- 
cal calculations of the average sign and comparison with 
the discrete-time QMC scheme are presented in the next 
Section. 



III. APPLICATIONS OF CT-QMC METHOD 

We test present algorithm for several well known mod- 
els in this Section. These examples show some of the 
advantages of the CT-QMC method. 

In all examples presented below we calculate a Green 
function at Matsubara frequencies G{iujn)- Total number 
of Matsubara frequencies is varied from 10 to 20. The 
typical number of QMC trials is 10^ 4- 10^. Normally, the 
error bar of the CT-QMC data for G{iujn) is less than 
3 • 10~^ for the lowest Matsubara frequency and becomes 
smaller as frequency increases. Obviously, values of these 
typical parameters depend on concrete system. 



A. Hubbard clusters 

To test the scheme, we start from a single isolated Hub- 
bard atom and a 2 x 2 Hubbard cluster. Results are 
compared with the known exact solution (see e.g. Ref. 
0). 

The solution for the atomic limit reads as follows: 

G(iLj) = -1^ + . " (35) 

Results for U = 2,/3 = 16,^ = U/2 are presented in Fig- 
ure|21 Thus CT-QMC data are in an excellent agreement 
with the analytical solution. 

Further we apply CT-QMC algorithm to the 2 x 2 Hub- 
bard lattice to compare with the auxiliary-field quan- 
tum Monte Carlo scheme'^. We start with the half- 
filled case — U/2, four electrons in the system). It 
can be shown that for the particular case of half-filling 
one can choose a^- — — 0.5 due to the particle- 
hole symmetry. Expression © for this case becomes 
< k >= f3N{0.5- < n-^ni >) with = 4. It can 
be verified that this choice delivers the minimal possi- 
ble < k >. Series © contains only the terms with an 
even k in this case, so it's appropriate to use steps ±2. 
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FIG. 2: Imaginary part of the Green function at Matsub- 
ara frequencies for a single atom with Hubbard repulsion U. 
Symbols are CT-QMC data, line is an exact solution^'. Pa- 
rameters: U — 2, P = 16,^1 — U/2. Error bar is less than 
symbol size. 
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FIG. 3: 2x2 Hubbard lattice at half-filling. Imaginary part 
of the Green function at Matsubara frequencies: symbols are 
CT-QMC data, line is exact-diagonalization data. Parame- 
ters: U = 4:,t — 1, l3 — 8, fj, = U/2. Error bar is less than 
symbol size. 



Results for C/ = 4, i = 1, /3 = 8 in comparison with the 
exact-diagonalization data are shown in Figure 13 

Cases of a single atom and a half-filled cluster do not 
suffer a sign problem. One can discuss a sign problem 
considering 2x2 Hubbard lattice away from half-filling. 
For this case a choice for a's was used. We concen- 
trate on the worst sign-problem case when there are three 
electrons in the systemSi. The average sign is presented 
in Figure ^ as a function of inverse temperature /3. We 
would like to stress that the CT-QMC algorithm agrees 
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FIG. 4: 2x2 Hubbard lattice away from half-filling: three 
electrons in the system. Average sign as a function of /3: 
CT-QMC (filled symbols) and auxiliary-field quantum Monte 
Carlc^ (opened symbols) algorithms results. Lines are guides 
to the eye. Parameters: U = 4, t = 1. 




FIG. 5: Real and imaginary parts of the Green function for 
2x2 Hubbard lattice away from half-filling: three electrons in 
the system. Parameters: U = 4, i = 1, /? = 14. Symbols are 
CT-QMC data, lines are exact-diagonalization results. Error 
bar for iu) > 2 is less than symbol size. 

with the auxiliary-field quantum Monte Carlo^ scheme 
(Figure El- Even for a relatively small average sign, nu- 
merical data remain to be in a good agreement with the 
exact-diagonalization, as Figure |5l shows. 



B. Metal-insulator transition on the Bethe lattice 

One of the advantages of the CT-QMC algorithm is a 
possibility to perform simulations at lower temperatures 



with higher accuracy than the auxiliary-field quantum 
Monte Carloi method. Here we present results for the 
metal-insulator phase transition in Hubbard model on 
Bethe lattice-. The effective one-site problem based on 
the dynamical mean-field theory'' is solved by CT-QMC 
method. 

The standard self-consistent loop of DMFT equations 
is as follows^. One starts with some initial guess for the 
Green function which is used to obtain the local Green 
function Q from the effective action a.sm 

g{r,T')=<Tcl,c^>s^„ig,). (36) 

A new guess for the Green function Qq is obtained from 
the equation for Bethe lattice {t — 1/2)1: 

Qf^^{iuj) — iuj + ^ - t'^Q{iuj). (37) 

Formulas (|36I37|I form a self-consistent loop of DMFT 
equations. The Green function which corresponds to the 
semi-circular density of states with band-width 2 is usu- 
ally used for the Bethe lattice: 

GoHoj) = pJ= . (38) 

The self-energy Yi{iLo) can be obtained from the follow- 
ing formula after the iteration procedure for the DMFT 
equations (|36I37|I has converged: 

Ys{iuj)^g~\iuj)-g-\ioj). (39) 

Results for the metal-insulator phase transition in 
Hubbard model on Bethe lattice at half-filling for /? = 64 
are presented in FigureEl Local Green functions and cor- 
responding self-energies are shown for values of Coulomb 
interaction U from the value [/ = 2 to the value U — 2> 
with the step AJ7 = 0.2. The results show a phase tran- 
sition from the metallic state (smaller values of U) to 
the insulating state (larger values of U) with a coexis- 
tence region in between. The data obtained agree well 
with previous studies of the transition where the stan- 
dard auxiliary-field quantum Monte Carloi algorithm 
was used as a solver for DMFT equations H36I37|> ^. Note, 
CT-QMC scheme gives better accuracy than auxiliary- 
field quantum Monte Carlch^- algorithm since one obtains 
the local Green function at Matsubara frequencies di- 
rectly in QMC. It allows one to perform simulations at 
lower temperatures. For instance, we tested the CT- 
QMC algorithm even at /? = 256 and obtained quite 
reasonable results for the metal-insulator phase transi- 
tion on Bethe lattice (see inset for Figure as well) . 

C. Multi-band model with a rotationally-invariant 
retarded exchange 

Another advantage of the CT-QMC algorithm is that 
it allows one to consider multi-band problems with inter- 
actions in the most general form: 

ijkl-.GG' 
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FIG. 6: Imaginary part of the Green function (a) and self- 
energy (b) at Matsubara frequencies for Hubbard model on 
Bethe lattice at half-filling obtained from the solution of self- 
consistent DMFT equations II36I37II by CT-QMC method. 
Parameters: /3 = 64, (7 = 2 3, AC/ = 0.2. AU data ob- 
tained with the initial guess for the Green function in the form 
(I38I I which corresponds to the metallic phase. Coexistence of 
metallic and insulating phases can be found, for example, at 
point U = 2.4. Inset shows data for the imaginary part of the 
Green function for /? = 256, U = 2.2 and U = 2.8. 



ant. The Gaussian part of the action represents the di- 
agonal semicircular density of states2^ with unitary half- 
band width H38(l . We used parameters C/ = 4, J = 1 at 
/3 = 4. Figures [7| and |S1 present the results for the lo- 
cal Green function G^s(r) and the four-point correlator 



,Oir„t 



iC 



ITr' 



>. The later quantity 



characterizes the spin-spin correlations and would vanish 
if the exchange were absent. 

A modification of this model was also studied where 
spin-flip operators were replaced with the terms fully 
non-local in time. For example, operator c^q-^^c^^'^ c\t^^c^'^'^ 



OT-r 



^Oir t „1Tt' 



As it is 



was replaced with i3~^ J dr'c, 

pointed in the introduction, the retardation effects in the 
interaction always appear if certain non-Gaussian degrees 
of freedom are integrated out. Therefore it is of impor- 
tance to demonstrate that CT-QMC scheme is able to 
handle the retarded interaction. 

The Green function in the time domain was obtained 
by a numerical Fourier-transform from the CT-QMC 
data for G{iuJn)- For high harmonics the following 
asymptotic form was used: — Im(zw + e)^^ with e « 2.9. 
The obtained dependencies are presented in Figure \7\ 
Results for the local and non-local in time spin-flip inter- 
actions are shown with solid and dot lines, respectively. 
It is interesting to note that the Green function is rather 
insensitive to the details of spin-flip retardation. The 
maximum-entropy guess for DOS is presented in the in- 
set to Figure [7| Both Green functions are very similar 
and correspond to qualitatively the same density of states 
(DOS). 

To demonstrate the effects due to retardation we cal- 
culated the four-point quantity x(t). These data are ob- 
tained similarly, the difference is that xi^'-^) is defined at 
Bose Matsubara frequencies and obeys a l/w^ decay. It 
turns out that a switch to the non-local in time exchange 
modifies x(t) dramatically. The local in time exchange 
results in a pronounced peak of x(t) at t « 0, whereas 
the non-local spin-fiip results in almost time- independent 
spin-spin correlations (Figure ISJ. 



We apply the proposed CT-QMC for the impor- 
tant problem of the super-symmetric two band impurity 
model at half-flUingSSiS^. To our knowledge, this is the 
first successful attempt to take the off-diagonal exchange 
terms of this model into account. These terms are im- 
portant for the realistic study of the multi-band Kondo 
problem because they are responsible for the local mo- 
ment formations^. The interaction in this model has the 
following form 

|(iV(r)-2)(7V(r)-2)-^(S(r).S(r)-fL(r).L(r)), (41) 

where N is the operator of total number, S and L are to- 
tal spin and orbital-momentum operators, respectively. 
The interaction is spin- and orbital- rotationally invari- 



IV. CONCLUDING REMARKS 

In conclusion, we have developed a fcrmionic continu- 
ous time quantum Monte Carlo method for general non- 
local in space-time interactions. It's successfully tested 
for a number of models. 

We demonstrated that for Hubbard-type models the 
computational time for a single trial step scales similarly 
to that for the schemes based on a Stratonovich transfor- 
mation. An important difference occurs however for the 
non-local interactions. Consider, for example, a system 
with a large Hubbard U and much smaller but still im- 
portant Coulomb interatomic interaction. One needs to 
introduce TV^ auxiliary fields per time slice instead of N 
to take the long-range forces into account. On the other 
hand, the complexity of the present algorithm remains 
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FIG. 7: Imaginary-time Green Function for the rotationally- 
invariant two-band model. Solid and dot lines correspond to 
the static and to the nonlocal in time spin-flip, respectively. 
The inset shows DOS estimated from the Green function. 




FIG. 8: Imaginary-time dependence of the four-point quan- 
tity xi^ ~ ''"') cJi^c^'^^cJj^^jC"'^^^ > for the rotationally- 
invariant two-band model. Solid and dot lines correspond to 
the static and to the nonlocal in time spin-flip, respectively. 



almo st the same as for the local interactions, because 
\W\ does not change much. This should be useful for the 
realistic cluster DMFT calculations and for the apphca- 
tions to quantum chemistry^. It is also possible to study 
the interactions retarded in time, particularly the super- 
exchange and the effects related to dissipation. This was 
demonstrated for an important case of the fully rotation- 
ally invariant multi-band model and its extension with 
non-local in time spin-flip terms. 



For the case of the Hubbard model the sign problem 
was found to be similar to what occurs for the auxiliary- 
field quantum Monte Carlo^ scheme. Nevertheless a gen- 
eral time-dependent form of the action (Eq.®) opens, in 
principal, the possibility for a two-stage renormalization 
treatment. Suppose we know a certain renormahzation 
of action, based on the local DMFT-solution as a starting 
point. Since DMFT is already a very good approxima- 
tion, we can expect the thus renormalized interaction to 
be smaller than the initial one, although it is perhaps 
nonlocal in time. Then one could expect that the lat- 
tice calculations with a renormalized interaction show a 
smaller sign problem. Practical investigation of such con- 
structed renormalization is a subject of the future work. 
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